From NYC Planning:

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.5     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.0.2     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.1
## Warning: package 'tibble' was built under R version 4.1.1
## Warning: package 'tidyr' was built under R version 4.1.1
## Warning: package 'readr' was built under R version 4.1.1
## Warning: package 'purrr' was built under R version 4.1.1
## Warning: package 'dplyr' was built under R version 4.1.1
## Warning: package 'stringr' was built under R version 4.1.1
## Warning: package 'forcats' was built under R version 4.1.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(tidycensus)
## Warning: package 'tidycensus' was built under R version 4.1.1
library(jsonlite)
## Warning: package 'jsonlite' was built under R version 4.1.1
## 
## Attaching package: 'jsonlite'
## The following object is masked from 'package:purrr':
## 
##     flatten
library(magrittr)
## Warning: package 'magrittr' was built under R version 4.1.1
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
library(sf)
## Warning: package 'sf' was built under R version 4.1.1
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
library(tigris)
## Warning: package 'tigris' was built under R version 4.1.1
## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
## 
## Attaching package: 'tigris'
## The following object is masked from 'package:tidycensus':
## 
##     fips_codes
library(censusapi)
## Warning: package 'censusapi' was built under R version 4.1.1
## 
## Attaching package: 'censusapi'
## The following object is masked from 'package:methods':
## 
##     getFunction
readRenviron("~/.Renviron")

get_decennial(
  geography = "tract",
  state = "NY",
  variables = "P1_001N",
  year = 2020,
  geometry = T
)
## Getting data from the 2020 decennial Census
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## Using the PL 94-171 Redistricting Data summary file
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================| 100%
## Note: 2020 decennial Census data use differential privacy, a technique that
## introduces errors into data to preserve respondent confidentiality.
## i Small counts should be interpreted with caution.
## i See https://www.census.gov/library/fact-sheets/2021/protecting-the-confidentiality-of-the-2020-census-redistricting-data.html for additional guidance.
## This message is displayed once per session.
#load_variables(2020, "pl")

varlist <- fromJSON("https://api.census.gov/data/2020/dec/pl/variables.json", simplifyDataFrame = T) %>% 
  extract2("variables") %>% 
  as.data.frame() %>% 
  mutate_all(~as.character(.)) %>% 
  pivot_longer(everything(), names_sep = "\\.",
               names_to = c("varname", "type"))

Load census data

ny_t_2020 <- st_read(dsn = "C:/Users/jg6849/Documents/Github/census-2020/data/tl_2020_36_tract/tl_2020_36_tract.shp") %>% 
  filter(COUNTYFP %in% c("005", "047", "061", "081", "085"))
## Reading layer `tl_2020_36_tract' from data source 
##   `C:\Users\jg6849\Documents\Github\census-2020\data\tl_2020_36_tract\tl_2020_36_tract.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5411 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.76259 ymin: 40.47658 xmax: -71.77749 ymax: 45.01586
## Geodetic CRS:  NAD83
ny_t_2020 
ny_t_2010 <- st_read(dsn = "C:/Users/jg6849/Documents/Github/census-2020/data/tl_2010_36_tract10/tl_2010_36_tract10.shp") %>% 
  filter(COUNTYFP10 %in% c("005", "047", "061", "081", "085")) %>% 
  rename(GEOID = GEOID10)
## Reading layer `tl_2010_36_tract10' from data source 
##   `C:\Users\jg6849\Documents\Github\census-2020\data\tl_2010_36_tract10\tl_2010_36_tract10.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 4919 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.76259 ymin: 40.4774 xmax: -71.77749 ymax: 45.01586
## Geodetic CRS:  NAD83
ny_t_2010
#ny_t_2020 <- tracts("NY", county = c("New York", "Kings", "Queens", "Richmond", "Bronx"),
#                year = 2020)
# 
# ny_t_2010 <- tracts("NY", county = c("New York", "Kings", "Queens", "Richmond", "Bronx"),
#                year = 2010)
# ny_t_2020
# ny_t_2010
ny_t_2020_filt <- ny_t_2020 %>%
  filter(AWATER < ALAND)

ny_t_2010_filt <- ny_t_2010 %>%
  filter(AWATER10 < ALAND10)

tract_intersection <- function(gdf1, gdf2, id) {
  gdf1 %>%
    filter(GEOID == id) %>%
    st_intersection(filter(gdf2, GEOID == id))
}

tract_difference <- function(gdf1, gdf2, id) {
  gdf1 %>%
    filter(GEOID == id) %>%
    st_difference(filter(gdf2, GEOID == id))
}

ny_t_int <- map_dfr(unique(ny_t_2020_filt$"GEOID"), ~tract_intersection(ny_t_2020_filt, ny_t_2010_filt, .))
ny_t_int
ny_t_diff_2020 <- map_dfr(unique(ny_t_2020_filt$"GEOID"), ~tract_difference(ny_t_2020_filt, ny_t_2010_filt, .))
ny_t_diff_2020
ny_t_diff_2010 <- map_dfr(unique(ny_t_2010_filt$"GEOID"), ~tract_difference(ny_t_2010_filt, ny_t_2020_filt, .))
ny_t_diff_2010
stopifnot(st_crs(ny_t_2020_filt) == st_crs(ny_t_2020_filt))

Intersection of census tracts between 2020 and 2010

ggplot(data = ny_t_int) +
  geom_sf(fill = "gray", lwd = 0.5) +
  theme_classic() +
  theme(
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(), 
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

Difference in census tracts between 2020 and 2010

ggplot() +
  geom_sf(data = ny_t_diff_2020, fill = "blue", color = "blue", alpha = 0.5) +
  geom_sf(data = ny_t_diff_2010, fill = "green", color = "green", alpha = 0.5) +
  theme_classic() +
  theme(
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(), 
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

ny_t_diff_comp <- ny_t_diff_2020 %>% 
  mutate(area_2020diff = st_area(geometry)) %>% 
  as.data.frame() %>% 
  select(c("GEOID", "area_2020diff")) %>% 
  full_join(ny_t_diff_2010 %>% 
              mutate(area_2010diff = st_area(geometry)) %>% 
              as.data.frame() %>% 
              select(c("GEOID", "area_2010diff")), by = "GEOID") %>% 
  left_join(ny_t_2020_filt, ., by = "GEOID") %>% 
  mutate(area_2020_total = as.numeric(st_area(geometry))) %>% 
  group_by(GEOID) %>% 
  mutate(area_ratio = as.numeric((area_2010diff + area_2020diff) / area_2020_total)) %>% 
  arrange(desc(area_ratio)) %>% 
  select("area_ratio", everything()) %>% 
  ungroup() %>% 
  mutate(area_ratio = case_when(is.na(area_ratio) ~ 0,
                                TRUE ~ area_ratio))

ny_t_diff_comp
stopifnot(eeptools::isid(ny_t_diff_comp, "GEOID"))

Scored difference

ggplot() +
  geom_sf(data = ny_t_diff_comp, aes(fill = area_ratio)) +
  theme_classic() +
  theme(
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(), 
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  scale_fill_continuous(high = "#132B43", low = "#56B1F7")

library(leaflet)
## Warning: package 'leaflet' was built under R version 4.1.1
library(mapview)
## Warning: package 'mapview' was built under R version 4.1.1
#m <- mapview(ny_t_diff_comp)
#m
ny_t_diff_comp <- ny_t_diff_comp %>% 
  filter(!is.na(area_ratio)) %>% 
  st_transform(crs = st_crs(4326)) %>% 
  filter(area_ratio > 0)

pal_rev <- colorNumeric(colorRamp(c("#40899A", "#FFFFFF"), 
                              interpolate="spline"), 
                    domain = ny_t_diff_comp$area_ratio)

pal <- colorNumeric(colorRamp(c("#40899A", "#FFFFFF"), 
                              interpolate="spline"), 
                    domain = ny_t_diff_comp$area_ratio, reverse = T)
pal
## function (x) 
## {
##     if (length(x) == 0 || all(is.na(x))) {
##         return(pf(x))
##     }
##     if (is.null(rng)) 
##         rng <- range(x, na.rm = TRUE)
##     rescaled <- scales::rescale(x, from = rng)
##     if (any(rescaled < 0 | rescaled > 1, na.rm = TRUE)) 
##         warning("Some values were outside the color scale and will be treated as NA")
##     if (reverse) {
##         rescaled <- 1 - rescaled
##     }
##     pf(rescaled)
## }
## <bytecode: 0x000000000c6281f8>
## <environment: 0x0000000004c11240>
## attr(,"colorType")
## [1] "numeric"
## attr(,"colorArgs")
## attr(,"colorArgs")$na.color
## [1] "#808080"
leaflet(data = ny_t_diff_comp) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addPolygons(data = ny_t_2010_filt,
              color = "black", stroke = T, weight = 0.2,
              fill = T, fillColor = "grey",
              group = "Census Tracts 2010") %>% 
  addPolygons(color = "gray", stroke = T,weight = 0.2,
              fill = T, fillColor = ~pal(area_ratio), fillOpacity = 1,
              group = "Census Tracts 2020") %>% 
  addLegend(position = "topright",
            pal = pal_rev,
            values = ~area_ratio,
            title = "Area ratio",
            labFormat = labelFormat(transform = function(area_ratio) 
              sort(area_ratio, decreasing = TRUE))) %>% 
  addLayersControl(overlayGroups = c("Census Tracts 2010", "Census Tracts 2020"))
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'